#include "ms_distribution.h"
#ifndef PI
#define UNDEF_PI
#define PI 3.1415926535897932384626433
#endif
#include "fe_vecmat.cxx"



/* ******************************************************************* */
RandomDistribution::RandomDistribution(){
	gsl_rng_env_setup();
	T = gsl_rng_default; 
	r = gsl_rng_alloc (T);
}

void RandomDistribution::setSeed(unsigned long int s)
{
    gsl_rng_set (r, s);
}

/* return a uniformly distributed scalar value on [0, 1] */
double RandomDistribution::RandVal() const{
	return gsl_rng_uniform(r);
	//hier zufallswert generieren
}

/* return a vector of uniformly distributed scalar value on [0, 1] */
void RandomDistribution::RandVal( double * x, const int N ) {
	for (int i=0;i<N;i++)	{
		x[i]=RandVal();
	}
	/* i.O., aber besser: Bibliothek aufrufen und direkt N Werte anfordern */
}

RandomDistribution::~RandomDistribution(){ 
	gsl_rng_free(r);
}


/* ******************************************************************* */
//UniformDistribution
UniformDistribution::UniformDistribution( const double aa, const double bb ){
	if( b > a )	{	a	= aa;	b	= bb; }
	else		{	b	= aa;	a	= bb; }
}

double UniformDistribution::GetRandomValue() const{
	return RandVal()*(b-a)+a;
}

void UniformDistribution::GetRandomVector(double * x, const int N ) const{
	/* besser: gleichverteilt N Werte 0..1, danach Schleife wie s.o. */
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
	}
}

/* ******************************************************************* */
//GaussianDistribution

GaussianDistribution::GaussianDistribution( const double _sigma ) {
	sigma	= _sigma; // Standard Deviation
}


double	GaussianDistribution::GetRandomValue() const 
{
	return gsl_ran_gaussian ( r ,  sigma);
	//Use the transformation z = mu + x on the numbers = ? //	
	//returned by gsl_ran_gaussian to obtain a Gaussian distribution with mean mu.//
}

void	GaussianDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
	}
}

/* ******************************************************************* */
//LogNormalDistribution
/* The lognormal distribution has the form 

   p(x) dx = 1/(x * sqrt(2 pi sigma^2)) exp(-(ln(x) - zeta)^2/2 sigma^2) dx

   for x > 0. Lognormal random numbers are the exponentials of
   gaussian random numbers */

LogNormalDistribution::LogNormalDistribution( const double _zeta, const double _sigma ) {
	zeta	= _zeta; // Location Parameter
	sigma	= _sigma; // Scale Parameter
}


double	LogNormalDistribution::GetRandomValue() const {
	return gsl_ran_lognormal ( r , zeta , sigma);
}
void	LogNormalDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
	}
}


/* ******************************************************************* */
//ExponentialDistribution
/* The exponential distribution has the form

   p(x) dx = exp(-x/mu) dx/mu

   for x = 0 ... +infinity */


ExponentialDistribution::ExponentialDistribution( const double _mu) {
	mu	= _mu; // mean value
}


double	ExponentialDistribution::GetRandomValue() const {
	return gsl_ran_exponential ( r , mu );
	
}
void	ExponentialDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}


/* ******************************************************************* */
//LaPlaceDistribution
/* The two-sided exponential probability distribution is  
   p(x) dx = (1/(2 a)) * exp(-|x/a|) dx
   for -infinity < x < infinity. It is also known as the Laplace distribution.  */

LaPlaceDistribution::LaPlaceDistribution( const double _a) {
	a	= _a; // scale parameter
}


double	LaPlaceDistribution::GetRandomValue() const {
	return gsl_ran_laplace ( r , a );
}
void	LaPlaceDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//CauchyDistribution
/* The Cauchy probability distribution is 

   p(x) dx = (1/(pi a)) (1 + (x/a)^2)^(-1) dx

   It is also known as the Lorentzian probability distribution */


CauchyDistribution::CauchyDistribution( const double _a) {
	a	= _a; // scale parameter
}


double	CauchyDistribution::GetRandomValue() const {
	return gsl_ran_cauchy ( r , a );
/*	in code:
gsl_ran_cauchy (const gsl_rng * r, const double a)
	{
		double u;
		do
		{
			u = gsl_rng_uniform (r);
		}
		while (u == 0.5);

	return a * tan (M_PI * u);
	}*/

	
}
void	CauchyDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//RayleighDistribution
/* The Rayleigh distribution has the form

   p(x) dx = (x / sigma^2) exp(-x^2/(2 sigma^2)) dx

   for x = 0 ... +infinity */


RayleighDistribution::RayleighDistribution( const double _sigma) {
	sigma	= _sigma; // scale parameter
}


double	RayleighDistribution::GetRandomValue() const {
	return gsl_ran_rayleigh ( r , sigma );
}

void	RayleighDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//LandauDistribution
/* The distribution is given by the complex path integral,
 *
 *  p(x) = (1/(2 pi i)) \int_{c-i\inf}^{c+i\inf} ds exp(s log(s) + x s) 
 *
 * which can be converted into a real integral over [0,+\inf]
 *
 *  p(x) = (1/pi) \int_0^\inf dt \exp(-t log(t) - x t) sin(pi t)
 */

LandauDistribution::LandauDistribution() {
	// no arguments required
}


double	LandauDistribution::GetRandomValue() const {
	return gsl_ran_landau ( r );

	
}
void	LandauDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
	}
}



/* ******************************************************************* */
//LevyAlphaStableDistribution
/* The stable Levy probability distributions have the form
 *
 * p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)
 *
 * with 0 < alpha <= 2. 
 *
 * For alpha = 1, we get the Cauchy distribution
 * For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c. */


LevyAlphaStableDistribution::LevyAlphaStableDistribution(const double _c, const double _alpha) {
	c	= _c; // scale
	alpha	= _alpha; // power (1=cauchy, 2=gaussian)

}


double	LevyAlphaStableDistribution::GetRandomValue() const {
	return gsl_ran_levy ( r , c , alpha);

	
}
void	LevyAlphaStableDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
	}
}




/* ******************************************************************* */
//LevySkewAlphaStableDistribution
/* The stable Levy probability distributions have the form
 * 
 * 2*pi* p(x) dx
 * 
 *   = \int dt exp(mu*i*t-|sigma*t|^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))) for alpha!=1
 *   = \int dt exp(mu*i*t-|sigma*t|^alpha*(1+i*beta*sign(t)*2/pi*log(|t|)))   for alpha==1
 * 
 * with 0<alpha<=2, -1<=beta<=1, sigma>0.
 * 
 * For beta=0, sigma=c, mu=0, we get gsl_ran_levy above.
 *
 * For alpha = 1, beta=0, we get the Lorentz distribution
 * For alpha = 2, beta=0, we get the Gaussian distribution

*/

 LevySkewAlphaStableDistribution::LevySkewAlphaStableDistribution(const double _c, const double _alpha, const double _beta) {
	c	= _c; // scale
	alpha	= _alpha; // power (1= cauchy, 2 = gaussian)
	beta    = _beta; // skew
}


double	LevySkewAlphaStableDistribution::GetRandomValue() const {
	return gsl_ran_levy_skew ( r , c , alpha, beta);
	
}
void	LevySkewAlphaStableDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
	}
}

/* ******************************************************************* */
//GammaDistribution
//  The Gamma distribution of order a>0 is defined by:
// 
//    p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx
// 
//    for x>0. 

GammaDistribution::GammaDistribution( const double aa, const double bb ){
	a	= aa; // a order
	b	= bb; // b scale
}

double GammaDistribution::GetRandomValue() const{
	return gsl_ran_gamma ( r , a , b );
}

void GammaDistribution::GetRandomVector(double * x, const int N ) const{
	/* besser: gleichverteilt N Werte 0..1, danach Schleife wie s.o. */
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//ChisquaredDistribution
/* The chisq distribution has the form

   p(x) dx = (1/(2*Gamma(nu/2))) (x/2)^(nu/2 - 1) exp(-x/2) dx

   for x = 0 ... +infinity */

ChiSquaredDistribution::ChiSquaredDistribution( const double _nu) {
	nu	= _nu; // degrees of freedom
}


double	ChiSquaredDistribution::GetRandomValue() const {
	return gsl_ran_chisq ( r , nu );

	
}
void	ChiSquaredDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//FDistribution
/* The F distribution has the form

   p(x) dx = (nu1^(nu1/2) nu2^(nu2/2) Gamma((nu1 + nu2)/2) /
   Gamma(nu1/2) Gamma(nu2/2)) *
   x^(nu1/2 - 1) (nu2 + nu1 * x)^(-nu1/2 -nu2/2) dx */

 FDistribution::FDistribution( const double _nu1 ,const double _nu2) { 
	nu1	= _nu1; //
	nu2	= _nu2; //
}


double	FDistribution::GetRandomValue() const {
	return gsl_ran_fdist ( r , nu1 , nu2 );

	
}
void	FDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//TDistribution
/* The t-distribution has the form

   p(x) dx = (Gamma((nu + 1)/2)/(sqrt(pi nu) Gamma(nu/2))
   * (1 + (x^2)/nu)^-((nu + 1)/2) dx */

TDistribution::TDistribution( const double _nu) {
	nu	= _nu; // degrees of freedom
}


double	TDistribution::GetRandomValue() const {
	return gsl_ran_tdist ( r , nu );

	
}
void	TDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}


/* ******************************************************************* */
//BetaDistribution
/* The beta distribution has the form

   p(x) dx = (Gamma(a + b)/(Gamma(a) Gamma(b))) x^(a-1) (1-x)^(b-1) dx */

BetaDistribution::BetaDistribution( const double aa, const double bb ){
	if( b > a )	{	a	= aa;	b	= bb; } // a, b: shape parameters
	else		{	b	= aa;	a	= bb; }
}

double BetaDistribution::GetRandomValue() const{
	return gsl_ran_beta ( r , a , b );
}

void BetaDistribution::GetRandomVector(double * x, const int N ) const{
	/* besser: gleichverteilt N Werte 0..1, danach Schleife wie s.o. */
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//LogisticDistribution
/* The logistic distribution has the form,

   p(x) dx = (1/a) exp(-x/a) / (1 + exp(-x/a))^2 dx

   for -infinity < x < infinity */


LogisticDistribution::LogisticDistribution( const double aa) {
	a	= aa;
}


double	LogisticDistribution::GetRandomValue() const {
	return gsl_ran_logistic ( r , a );

	
}
void	LogisticDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}

/* ******************************************************************* */
//WeibullDistribution
/* The Weibull distribution has the form,

   p(x) dx = (b/a) (x/a)^(b-1) exp(-(x/a)^b) dx

 */

WeibullDistribution::WeibullDistribution( const double aa, const double bb ){
	a = aa; // scale parameter
	b = bb; // exponent
}

double WeibullDistribution::GetRandomValue() const{
	return gsl_ran_weibull ( r , a , b);
}

void WeibullDistribution::GetRandomVector(double * x, const int N ) const{
	/* besser: gleichverteilt N Werte 0..1, danach Schleife wie s.o. */
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
	}
}



/* ******************************************************************* */
//Gumbel1Distribution
/* The Type I Gumbel distribution has the form,

   p(x) dx = a b exp(-(b exp(-ax) + ax)) dx */

Gumbel1Distribution::Gumbel1Distribution( const double aa, const double bb ){
	a = aa; // order
	b = bb; // scale parameter
}

double Gumbel1Distribution::GetRandomValue() const{
	return gsl_ran_gumbel1 ( r , a , b);
}

void Gumbel1Distribution::GetRandomVector(double * x, const int N ) const{
	/* besser: gleichverteilt N Werte 0..1, danach Schleife wie s.o. */
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}


/* ******************************************************************* */
//Gumbel2Distribution
/* The Type II Gumbel distribution has the form,

   p(x) dx = b a x^-(a+1) exp(-b x^-a)) dx

 */

Gumbel2Distribution::Gumbel2Distribution( const double aa, const double bb ){
	a = aa; // order
	b = bb; // scale parameter
}

double Gumbel2Distribution::GetRandomValue() const{
	return gsl_ran_gumbel2 ( r , a , b );
}

void Gumbel2Distribution::GetRandomVector(double * x, const int N ) const{
	/* besser: gleichverteilt N Werte 0..1, danach Schleife wie s.o. */
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}


/* ******************************************************************* */
//PoissonDistribution
/* The poisson distribution has the form

   p(n) = (mu^n / n!) exp(-mu) 

   for n = 0, 1, 2, ... .*/

PoissonDistribution::PoissonDistribution( const double _mu ) {
	mu	= _mu; // scale parameter
}


double	PoissonDistribution::GetRandomValue() const {
	return gsl_ran_poisson ( r , mu );
}

void	PoissonDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}



/* ******************************************************************* */
//BinomialDistribution
/* The binomial distribution has the form,

   prob(k) =  n!/(k!(n-k)!) *  p^k (1-p)^(n-k) for k = 0, 1, ..., n */


BinomialDistribution::BinomialDistribution( const double _p , const unsigned int _n ) {
	p	= _p; // probability
	n	= _n; // Number of trials
}


double	BinomialDistribution::GetRandomValue() const {
	return gsl_ran_binomial ( r , p , n );
}

void	BinomialDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}


/* ******************************************************************* */
//PascalDistribution
/* The Pascal distribution is a negative binomial with valued integer n

   prob(k) =  (n - 1 + k)!/(n!(k - 1)!) *  p^n (1-p)^k for k = 0, 1, ..., n

   */

PascalDistribution::PascalDistribution( const double _p , const unsigned int _n ) {
	p	= _p; // probability
	n	= _n; // order (integer)
}


double	PascalDistribution::GetRandomValue() const {
	return gsl_ran_pascal ( r , p , n );
}

void	PascalDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}


/* ******************************************************************* */
//LogarithmicDistribution
/* Logarithmic distribution 

   prob(n) =   p^n / (n log(1/(1-p)) for n = 1, 2, 3, ... */

LogarithmicDistribution::LogarithmicDistribution( const double _p) {
	p	= _p; //probability
}


double	LogarithmicDistribution::GetRandomValue() const {
	return gsl_ran_logarithmic ( r , p );
}

void	LogarithmicDistribution::GetRandomVector(double * x, const int N ) const  {
	for (int i=0;i<N;i++) {
		x[i]=GetRandomValue();
		}
}


/* ******************************************************************* */


MisesFisherDistribution::MisesFisherDistribution(){

}



/* ******************************************************************* */


Histogram::Histogram( const Vector & x, const Vector & W, const int N, const bool store ) {
	has_dat=store;
	dat		= x;
	n_bin	= N;
	Recalculate();
	gsl_h = 0;
	if (store==false) {
		/* compute only the histogram, but do not copy x to dat */
		dat.~Vector(); /* call destructor of vector class */
	}
}

/* Gewichtung ist in diesem Fall ein Vektor... */
/* copy vector W, since values are changed!!! */


Vector Histogram::GetW(Vector W) {
// 	double variance = gsl_stats_variance(dat, 1, n_bin);
	double SummeWW = 0.; /* initialize value */
	for (int i=0;i<n_bin;i++){
		SummeWW += W[i]; 
	}
	
	W = W / SummeWW; // entspricht W = W / SummeWW
	// W.X; /* double * auf die Daten von W */
	return W;
};


Histogram::~Histogram() {
	gsl_histogram_free(gsl_h);
}

void Histogram::Recalculate() {
	hist.resize(n_bin,true);
	if (dat.n<1)
		return;
	v_min=dat[0];
	v_max=dat[0];
	for (int i=1;i<dat.n;i++) {
		if (dat[i]<v_min)
			v_min=dat[i];
			//alternativ: Kapitel 21.8 der GSL: gsl_stats_min ( data , stride , n )
		else if (dat[i]>v_max)
			v_max=dat[i];
			//alternativ: Kapitel 21.8 der GSL: gsl_stats_max ( data , stride , n )
	}
	
	if( gsl_h != 0 ) gsl_histogram_free( gsl_h );
	h_bin	= (v_max-v_min)/(double)n_bin;
	
	gsl_h = gsl_histogram_alloc( n_bin );
	double range[ n_bin+1 ]; /* does not work for n_bin > several thousand */

		/* alternative implementation:
	for( int i=0; i<=n_bin; i++ ){
		range[i] = a + h_bin*(double)i;
	}
	gsl_histogram_set_ranges( gsl_h, range, n_bin+1 );	  */

	gsl_histogram_set_ranges_uniform( gsl_h, v_min, v_max );

	GetW(W);
	for (int j=0;j<dat.n;j++) {
		int interval=0; 
		gsl_histogram_accumulate( gsl_h, j, W.X[j]); 
		hist[interval]++;
	}
	
	WriteDatafile("MeinDateiName.txt");
	return;
}

Vector	Histogram::GetHistogram() const {
	Vector h(hist);
	return h;
}

/*...*/


void	Histogram::SetHistogram( Vector &h) const {
	hist *h;
}


/* if has_dat --> clear dat and unset has_dat */
void	Histogram::ClearData() {
	if(has_dat) {
		dat = 0.;
		has_dat = false;
	}
}


/* set new data vector and recompute the histogram */
void	Histogram::SetData( const Vector & x, const bool store ) {
	dat *x;
	has_dat = true;
}

/* Return a copy of the data vector (if existing) into D and return true; otherwise: return false and delete D */
bool	Histogram::GetData( Vector & x ) const {
	return true;
}

/* set new number of bins and recompute; if has_dat = false --> return false (i.e. unable to change number of bins; keep old N) */
bool Histogram::SetBinNumber( const int _N ) {
	if (has_dat == false) {
		return false;
	} 
	else {
		int n_bin = _N;
	return true;
	}
}

// write Data File:
void  	Histogram::WriteDatafile( const char fn[] ) const {
char acFileName[] = "MeinDateiName.txt";
WriteDatafile(acFileName);
//alternativ:
//WriteDatafile("MeinDateiName.txt");
//alternativ:
//char acFileName[] = "MeinDateiName.txt";
//char* pFileName = &(acFileName[4]); // pointer auf das "D"
//WriteDatafile(pFileName); // speichert dann in "DateiName.txt"
}


/* Write default GnuPlot script for the creation of EPS, PNG, ... figures */
void	 Histogram::WriteGnuplotScript( const char fn[] )  {
	ofstream scriptfile;
	scriptfile.open(fn);
	scriptfile << "#\n";
	scriptfile << "#\n Head comment";
	scriptfile << "#\n";
	scriptfile << "set title \"default_title\"";
	
	/* etc. etc. */
	scriptfile.close();
}

/* Write default Matlab script for the creation of EPS, PNG, ... figures */	
void	 Histogram::WriteMatlabScript( const char fn[] )  {
	ofstream scriptfile("Daten.mat"); 
	for(int i = 0; i < dat.n; ++i){
		scriptfile << dat[i] << " " << dat[i+1] << "\n";
	}
	scriptfile.open(fn);
	scriptfile << "#\n";
	scriptfile << "#\n Head comment";
	scriptfile << "#\n";
	scriptfile << "set title \"default_title\"";

	/* etc. etc. */
	scriptfile.close(); 
	
}

/* ******************************************************************* */


